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Earthquakes are among the most destructive natural events. The 6 April 2009, 6.3-M„ earthquake in 
L'Aquila (Italy) markedly altered the karstic Gran Sasso Aquifer (GSA) hydrogeology and geochemistry. 
The GSA groundwater invertebrate community is mainly comprised of smaU-bodied, colourless, blind 
microcrustaceans. We compared abiotic and biotic data from two pre-earthquake and one post-earthquake 
complete but non-contiguous hydrological years to investigate the effects of the 2009 earthquake on the 
dominant copepod component of the obligate groundwater fauna. Our results suggest that the massive 
earthquake-induced aquifer strain biotriggered a flushing of groundwater fauna, with a dramatic decrease in 
subterranean species abundance. Population turnover rates appeared to have crashed, no longer 
replenishing the long-standing communities from aquifer fractures, and the aquifer became almost totally 
deprived of animal Ufe. Groundwater communities are notorious for their low resilience. Therefore, any 
major disturbance that negatively impacts survival or reproduction may lead to local extinction of species, 
most of them being the only survivors of phylogenetic lineages extinct at the Earth surface. Given the 
ecological key role played by the subterranean fauna as decomposers of organic matter and "ecosystem 
engineers", we urge more detailed, long-term studies on the effect of major disturbances to groundwater 
ecosystems. 
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Though groundwater is predominantly viewed as a life-devoid water reservoir, there is a unique, stable animal 
community which lives in subterranean waters' ^\ Life at extremes? Not for groundwater animals, which 
have developed adaptive traits as long-term groundwater dwellers'\ They are small-sized, colourless, blind, 
with low metabolic rate and fecundity, features which allow survival in total darkness, lack of photosynthesis, low 
food avaUabUity, and strict dependence on the organic matter entering groundwater from the surface"*'^. 
Crustaceans dominate in groundwater communities^, and most of them are minute copepods*" (Fig. la), followed 
by ostracods^ and larger-bodied amphipods*. 

More than 90% of the groundwater fauna is composed of local endemics, groundwater-exclusive taxa (from 
species through class) that significantly expand the phylogenetic global animal biodiversity""''"''''''. Despite its 
uniqueness, and a timid increase in presence within the scientific literature on inland waters, the groundwater 
fauna remains poorly investigated and the least known component in global biodiversity assessments"''" being 
hidden, invisible, and often unapproachable. Several groundwater taxa are living fossils, the only remnants of very 
ancient evolutionary lineages disappeared from the Earth surface, but preserved underground and at high risk of 
extinction^'"'". Groundwater fauna importance in the overall freshwater biodiversity goes beyond its uniqueness. 
Obligate-groundwater species (i.e. stygobites, which complete their life cycle in groundwater, being unable to 
reproduce and survive in surface-water habitats) play key roles in organic matter recycling'"''" and as "ecosystem 
engineers" because the burrowing life style of worm-like elements allows the oxygenation of the minute voids 
among sediment particles'*. Groundwater invertebrates may also serve as early warning systems in monitoring 
groundwater quality and quantity, being sensitive to groundwater contamination'^" and withdrawal'"", because 
they respond to such disturbances as much as surface-water species do. However, a knowledge vacuum stiU 
surrounds the response of groundwater fauna to both human-induced"" and natural drivers, earthquakes holding 
the primacy among the latter. 

Earthquakes are among the most destructive natural events'"". Their effects on underground hydrology are well 
known, and are associated with the potential loss of groundwater supply for human use^'. Hydrological responses 
depend on several factors, such as earthquake magnitude, distance from the epicentre, and local geological 
conditions"""^''. The mainshock of the L'Aquila earthquake on 6 April 2009 (M„ 6.3) dramatically changed the 
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— 1997 : mean annua! discharge (±SD): 

5.68 ± 0.21 m^s"' 

— 2005 : mean annua! discharge (±SD): 

6.02 ± 0.26 m^s"' 



2009 : mean annual discharge (±SD): 
5.85 ± 0.72 m^s'^ 



2012 : mean annual discharge (±SD): 
7.14 ±0.41 m^s'^ 



Figure 1 | Diagram of the response of groundwater biodiversity to changes of the Gran Sasso aquifer (central Italy) hydrodynamics before and after the 
earthquakeof 6 April 2009. SEM micrographs ofrepresentative groundwater crustacean copepods (white bar = 200 [im) (a), schematic representation of 
pre-seismic (b, e), co-seismic (c, f), and post-seismic (d, g) groundwater flow from surface recharge areas through the aquifer (main conduit, lateral 
chambers and fracture network) to the spring discharge area, with obligate groundwater microcrustacean responses. Bottom panels report the discharge 
of the Tirino Springs (TS) during the March-June periods for each sampling year; mean annual discharges ± standard deviations also are reported. 



groundwater flow of the Gran Sasso Aquifer (GSA)^"""^' and conse- 
quently of the Tirino River vaUey, where most of total aquifer dis- 
charge is located (Fig. 2). A strong base-flow component and the 
presence of no-flow limits defined by aquitards (constant head con- 
ditions) determine a steady water table at the Tirino Springs (TS), 
historically characterised by limited seasonal variations in spring- 
time recharge, summer peak and autumn/winter shortage^^"^^. 
Short-term earthquake effects on GSA hydrology were determined 
by a pore pressure increase due to aquifer deformation^^. This short- 
term process triggered a medium-term response in the fractured 
aquifer: an earthquake-induced increase in bulk hydraulic conduc- 
tivity in the recharge area, near the earthquake-ruptured fault zone, 
attributable to fracture clearing and/or microcrack formations'^, 
contributed to the anomalous rising of the water table and flow rate 
in discharge zones (Fig. 2). The otherwise stable water table at TS had 
progressively risen since the post-mainshock week to up to one metre 
(equivalent to lOX the historical data) in a month. Such an increase 
in water table levels could not be attributed to rainfall or snow- 
j^gj^25,27,29 -pj^g water table in December 2012 was stiU at least 1 m 
higher than in the two years after the mainshock. A partial return 
toward pre-earthquake discharge levels started only in 2013 
(Supplementary Fig. 1). The increased post-seismic discharge at TS 
was evident already in May 2009, reaching out-of-range values 
(>30% of the previous 15 years) until the end of 2012. 

A direct response of subterranean biodiversity to earthquakes 
could be expected, due to the tight relationship between ground 
shaking, aquifer strain, fracturing and fracture clearing and habitat 
changes for groundwater animals, but no earthquake-related bio- 



logical information in groundwater is available so far. What can 
groundwater fauna witness about the effects of earthquakes? How 
do they react to catastrophic events that dramatically rearrange the 
habitat where they spend their entire life cycle? Are they sensitive to 
earthquake-triggered hydrologic changes in groundwater? Where no 
other entrance to groundwater is available, answering these ques- 
tions requires entering into the aquifers through "open win- 
dows"'"'", such as springs'' or other natural (caves) or artificial 
underground openings (boreholes)'^. Springs may serve as a window 
to the underground, as they are considered "readily accessible out- 
crops of the planet's extensive groundwater resources"". 

We monitored microcrustacean copepod (body size at the |im 
scale) species richness and abundance along with physicochemical 
features at the TS, one of the main GSA discharge sites, before and 
after the earthquake to assess whether groundwater communities 
responded to hydrological changes triggered by the earthquake. 
Copepods were used as the target biological variable being by far 
the most abundant and species-rich groundwater group (—80% by 
individual-based density). The remainder of the meiofauna'' col- 
lected at TS included nematods (—14%), oligochaetes (—3%), ostra- 
cods (<1%), halacarid mites (<1%), gastropods (<1%), together 
with a small fraction of macroinvertebrates, such as amphipods 
and chironomid larvae (—3%). 

Results 

Changes in groundwater environmental descriptors. Univariate 
nonparametric tests returned significant 1997-2005-V5.-2012 diffe- 
rences for most of the environmental parameters analysed. Post- 
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Figure 2 | Study area. Location of the Gran Sasso Aquifer (GSA), its main recharge and discharge areas, including the Tirino Springs (TS) where sampling 
was carried out, and the epicentre of the earthquake of 6 April 2009 that hit L'Aquila and surroundings. The two red dots represent the main springs 
(Presciano and Capodacqua) within the TS system (software program CorelDRAW 11 (version 11.633), 2002 Corel Corporation). 



seismic groundwater pH was significantly lower and Ca^* and 
dissolved oxygen (DO) concentrations were significantly higher 
than pre-seismic levels (Fig. 3). Electrical conductivity increased in 
September 2009 (up to 580 |iS cm"')^^, but returned rapidly to pre- 
earthquake values (Fig. 3). Very fine sand and particulate organic 
matter (POM) significantly increased in 2012 at the TS outlets 
(Fig. 3). 

Water temperature at the TS sampling sites remained roughly 
constant throughout the sampling period, with median values at 
around 11°C (Supplementary Fig. 2). Post-seismic nitrate concentra- 
tions were significantly lower than pre-seismic values (Supple- 
mentary Fig. 2). Sulphate concentrations in 2012 returned to 1997 
values after a transient peak in 2005 (Supplementary Fig. 2). Except 
for very fine sand (Fig. 3), none of the grain size classes sampled at TS 
changed significantly across sampling years, though the finest classes 
(medium and fine sand, and sUt and clay) tended to increase after the 
earthquake, albeit qualitatively (Supplementary Fig. 2). The raw abi- 
otic data set is in Supplementary Table 1. 

Effects on groundwater biodiversity. Total groundwater copepod 
abundance declined significantly after the earthquake (Fig. 3). 
Though statistical separation remained incomplete, groundwater 
species richness also exhibited a post-seismic negative trend 
(Fig. 3). Community patterns for non-obligate groundwater 
copepods (i.e., species which live in spring habitats because cold 
stenothermic or generalists, and species originally coming from 
surface water in recharge areas and "leaked" underground) did not 
change significantly after the earthquake (Fig. 3). Though both non- 
obligate and obligate groundwater species were found at 26 out of the 



32 sampling sites in 2012, such numbers represented a decrease for 
the stygofauna and an increase for non-obligate groundwater 
copepods (Table 1). The response at species level was varied for 
surface-water taxa, with six (—46%) species decreasing and five 
(—38%) species increasing in site-specific occurrence in 2012, 
while seven of the nine obligate groundwater species (including all 
the three most common pre-earthquake species) had decreased 
significantly in post-earthquake site-specific occurrence (Table 1). 

Abundances of the 22 copepod species that we found at the 32 
sampling sites declined dramatically after the earthquake (Supple- 
mentary Table 2). Pre- vs. post-earthquake differences in species 
abundance were highly significant (Supplementary Table 3; pairwise 
tests after significant one-way PERMANOVA: 1997 vs. 2012: P = 
0.0009; 2005 vs. 2012: P = 0.0015), while differences between the two 
pre-seismic samplings remained non-significant (1997 vs. 2005: P = 
0.5131). Obligate groundwater species mostly accounted for the dra- 
matic post-earthquake decrease in faunal abundance (pairwise tests 
after significant one-way PERMANOVA: 1997 vs. 2012: P = 0.0007; 
2005 vs. 2012: P = 0.0012; 1997 vs. 2005: P = 0.7450), with the 
contribution from non-obligate groundwater species remaining 
non-significant (one-way PERMANOVA: P = 0.1379) (Supple- 
mentary Table 3). 

The groundwater species Nitocrella pescei, Diacyclops paolae, 
Parastenocaris lorenzae and Elaphoidella mabelae cumulatively pro- 
vided 90% of the discrimination between pre- vs. post-earthquake 
assemblages (Supplementary Table 4). N. pescei contributed the most 
(41.7%) to the pre- vs. -postdissimUaritybydecliningby — 2.5X after 
the earthquake (Fig. 4, Supplementary Table 4). The other three most 
common species contributed further to pre- vs. post-seismic dissim- 
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Figure 3 | Selected abiotic and biological variables recorded at the Tirino Springs before (1997, 2005) and after the 2009 earthquake (2012). 

Boxplot representations of median (thick black line), upper and lower quartiles (box), 1.5X interquartile range (whiskers) and outliers (circles) for 
selected abiotic (green/blue) and biological (yellow/orange) variables at TS in the three sampling years. Post-seismic values in blue or orange. Different 
letters within each panel represent significantly different data distributions [Nemenyi comparisons after significant Friedman tests (Q,P\df= 2 for 
each)]. 



ilarity to a cumulative —50%, though D. paolae and P. lorenzae 
decreased in abundance by only 0.5 X and E. mabelae had slightly 
increased by 2012 (Table 1, Supplementary Table 4, Fig. 4). The vast 
majority of the remaining groundwater species became rare to very 
rare in 2012 (Table 1, Supplementary Table 2). 

Discussion 

The mainshock of the L'AquUa earthquake of 6 April 2009 markedly 
changed the Gran Sasso groundwater flow and consequently the 
discharge at the main spring sites. Medium-term local effects 
included an anomalous increase in TS discharge through 2012"'^"^' 
(Supplementary Fig. 1). The TS discharge gradually returned toward 
pre-seismic values only in summer 2013 (Supplementary Fig. 1), 
consistent with model predictions^^'^'^. 

Pre- vs. post-earthquake changes in groundwater physico-chem- 
istry were observed at the TS main outlets. Post-seismic water pH 
(Fig. 3) decreased as a reflection of an increase of deeply-derived CO2 
along the main faults during the mainshock'''. The post-seismic 
decrease in pH was accompanied by a significant increase in Ca^* 
concentration from calcite dissolution^' (Fig. 3). Electrical conduc- 
tivity increased in September 2009 as a result of the clearing of 
previously obstructed fractures from which high residence-time 
ion-enriched water flushed out after crustal stress and ground shak- 
ing (Fig. 3). Concentrations of very fine sand and POM increased at 
the main TS outlets (Fig. Ib-d, Fig. 3) due to the massive dewatering 
of the conductive systems which emptied the annexed fractures in the 
aquifer. Dissolved oxygen (DO) concentration increased signifi- 
cantly after the earthquake (Fig. 3), as expected for the increase in 



fast groundwater flowpaths in the aquifer due to a post-seismic 
higher hydraulic conductivity. 

The decrease in post-seismic nitrate concentrations (Supplement- 
ary Fig. 2) was most likely due to a markedly lower human presence 
in the catchment area following the post-earthquake temporary 
evacuation of towns and farms. Higher sulphate concentrations in 
2005 (Supplementary Fig. 2) were related to the dual groundwater 
flow system feeding the TS''. The seasonal contribution from local 
heavy rainfall led to an above-average recharge rate in 2005, and to an 
enrichment in S04^" due to a more intense interaction with the 
Triassic evaporite bedrocks outcropping along the local and most 
superficial groundwater flowpath". Conversely, the highest post- 
seismic groundwater contribution was from the deepest, base-flow 
fractured aquifer, described by large karstic conduits and the annex 
capacitive subsystem^'"^". Absence of temporal changes for most 
such variables and the overall low values for nitrates and sulphates 
(Supplementary Fig. 2) suggest a limited (if any) influence of such 
variables on the observed changes in aquifer-dwelling biota. 

Stygobiotic copepods responded negatively to changes in the 
groundwater hydrodynamics triggered by the earthquake (Table 1, 
Figs. 3, 4, Supplementary Table 3). However, the relative contribu- 
tions of the single-event massive mainshock (from which the stygo- 
fauna had yet to recover as of late 2012) or from the post-mainshock 
diminished drift (due to the overall initial crash in stygobiotic abund- 
ance) remain unknown. Interestingly enough, our earthquake- 
targeting investigation has led to a noteworthy corollary finding: 
though most species within the harpacticoid genus Elaphoidella are 
true stygobites"", the slight but significant post-earthquake increase 
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Table 1 | Changes in frequency of occurrence of copepods in the Tirino Springs (TS). Number of sites per sampling year at which the two 
major groups of obligate groundwater (GW) and non-obligate groundwater (nGW) species, and all individual species, were recorded as 
present. Significant departures from expected values in 201 2 (according to Pearson adjusted residuals following significant contingency- 
table, Pearson-distance-based tests of independence on coded data: groups: Xp^ = 9.627, df= 2, P = 0.008; species: Xp^ = 5507.5 1 1 , 
c/f = 42, P< 0.0001 ) are highlighted with up and down arrows, representing significant increases or decreases, respectively; * = most 
common pre-seismic species 
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in E. mabelae presence (Table 1, Fig. 4) supports the view of an 
intermediate stygo- and crenobiotic ecology of this species^'. 

Though a few species can cope with the challenging fast-flowing 
conduit system of karstic aquifers, the dendritic fracture network 
preferentially allows the establishment of permanent (in ecological 
times) invertebrate communities'"*'^'^. The low flow rate in such frac- 
ture "chambers" (Fig. Ib-d) allows the storage of fine sediments. 



Nitocrella pescei 



Diacyclops paolae 




Figure 4 | Pre- and post-seismic abundance patterns of dominant 
groundwater copepods at tlie Tirino Springs. Nonmetric 
multidimensional scaling ordination of pre- and post- earthquake 
abundance of the obligate groundwater copepod species that accounted for 
—90% of abundance changes; 2D stress = 0.16. Pre-earthquake replicate- 
specific abundance data in light blue, and post-earthquake abundances in 
red; circle size is proportional to abundance). 



which form the physical substrate for burrowing animals, and of 
POM, which constitutes the bulk of the food. Open-water swimmers 
also avoid massive drifting (Fig. lb). 

Based on our limited dataset of three non-contiguous sampling 
years, and considering the relatively stable pre-seismic groundwater 
faunal patterns, we propose that the earthquake flushed out ground- 
water invertebrates, very fine sand and POM at the aquifer discharge 
sites (Fig. Ic). The massive earthquake-induced increase in hydraulic 
conductivity and the consequent aquifer dewatering (Fig. le-g) thus 
acted as what we call a biotrigger. Immediately after the mainshock, 
groundwater communities were stranded in springbed sediments, 
where they survived and maintained the pre-earthquake species rich- 
ness albeit at lower abundances, a condition expected if individuals of 
fracture-dwelling species were massively flushed out (Table 1). 
Neither the post-seismic higher DO nor the higher POM availability 
at the TS outlets favoured the survival and especially the reproduc- 
tion of groundwater species, and the contribution of post-earthquake 
drift for groundwater species was still negligible as of 2012. Stygobite 
numerical losses in the aquifer fractures were not consistently or fully 
counterbalanced by conspecifics coming from their original habitat 
(i.e. small and large fractures eventually less affected by the aquifer 
dewatering). Progressively decreasing abundances of obligate 
groundwater species (i.e., stygobites) were found at the TS outlet, 
because the aquifer became almost totally deprived of animal life, no 
longer replenishing groundwater populations normally flushing out 
at the spring outlets (Fig. Ic, d). Abundance and species richness of 
non-obligate groundwater species remained unchanged (Fig. 3), sug- 
gesting that the earthquake had an impact only on aquifer-dwelling 
microinvertebrates. We therefore believe that the main driver of the 
observed crash in post-earthquake biodiversity occurred inside the 
aquifer. 
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How much time is needed for the recolonisation of the Gran Sasso 
aquifer by groundwater species, allowing the reestablishment of 
viable microcrustacean populations? We think that such a process 
will take a very long time, probably several decades, much longer 
than the return to the pre-earthquake aquifer recharge/discharge 
cycle and geochemistry. Groundwater communities are notorious 
for their low resilience''. Virtually all species have low reproductive 
rates, long life spans (up to 3.5 years for groundwater copepods vs. 
£1 year for surface-water relatives), and long post-embryonic 
development times before reaching sexual maturity'* ''. The potential 
for replacement by other groundwater species is very low, aquifer 
communities being characterised by low community redundancy' 
(although omnivory is claimed as a relatively frequent functional 
trait in the oligotrophic groundwater''", an ecological statement 
awaiting confirmation). 

Our investigation is based on a limited dataset of three complete 
but non-contiguous hydrological years, so that a precise quantifica- 
tion of background variability ("noise") was not possible. Though 
some abiotic variables were fully comparable in the two pre- 
earthquake years (e.g., pH, Ca"'*, POM and DO in Fig. 3), others 
were not (electrical conductivity in Fig. 3 and SO/" in Supple- 
mentary Fig. 2). Also, some site-specific spatial variability occurred 
(e.g., electrical conductivity slightly higher at site 8 in Supplementary 
Table 1), suggesting the existence of some background variability in 
hydrogeology. However, all statistical analyses pointed to a pre-vs.- 
post-earthquake variation higher than within pre-earthquake vari- 
ation (e.g., Friedman tests in Fig. 3 and all multivariate biological 
tests). Therefore, we believe that our results remain valid and robust 
as a first-look investigation into the effects of a major single-event 
disturbance such as an earthquake on groundwater ecosystems, but 
we urge more detailed, long-term studies on the topic. Despite the 
limited dataset, our results strongly suggest that a catastrophic nat- 
ural disturbance such a major earthquake seems to negatively impact 
the survival or reproduction of aquifer- dwelling stygobites, with a 
potential detrimental effect on groundwater biodiversity. 

Methods 

Sampling and sample processing. The Tirino Springs (TS) were selected as the study 
site because a natural entrance to the aquifer (e.g., caves) is not available in the Gran 
Sasso karst. In such cases, sampling is carried out at spring outlets ^''^^'^^ or through 
artificial boreholes^^. Most of the Gran Sasso Aquifer (GSA) discharge occurs through 
springs at the aquifer boundary (Fig. 2), and the TS receive —65% of the aquifer 
discharge*^. We sampled the TS through bedrock fractures in the springbed, deep into 
the fractures well below the springbed itself, or in the aquitards where present, at 
different depths in order to collect obligate groundwater dwellers (i.e., stygobites). 
The high groundwater discharge from under the springbed maximised the collection 
of groundwater and stygobites, either drifted from the karstic discontinuities of the 
carbonate bedrock, or stranded in the alluvial aquitards. Stygobites depend on 
groundwater to complete their life cycle, as they cannot reproduce in other 
environments. Species were characterised as stygobites on the basis of well- 
established morphological, physiological and ecological features^-^-^'^^. Stygobites are 
frequently found in spring habitats where they can temporarily survive after being 
drifted from their original groundwater habitats^^'^^'*\ Surface-water species are also 
found along with stygobites, either flushed out from the aquifer where they enter from 
recharge areas, or because they live in springs, as crenobionts or surface-water habitat 
generalists^^. 

Abiotic and biological variables were determined at eight sampling sites at the TS 
adopting a random sampling method with four seasonal replicates per site (Jan-Feb, 
Mar-Apr, Jul-Aug, Sep-Oct) in three discrete years: 1997 and 2005 (before the 
earthquake) and 2012 (after the earthquake). Samples were collected with a Bou- 
Rouch pump*^ and mobile pipes hammered at each sampling site. Samples were 
composite from 0.3, 0.7, and 1.5 m below the springbed. Since both sample size and 
pumping rate may affect estimates of faunal density for a range of different sediment 
types^^'*^"^^, we maintained consistent pumping rates and sample volume sizes 
throughout the investigation. Though 5-L samples would suffice to detect spatial 
variation in invertebrate composition among sites with different sediment texture, 
sample size affecting mainly large-sized invertebrates like insect larvae*^ a sample size 
of 20 L was adopted for density estimates, thus including abundance-related rarity*^ 
of groundwater copepods*-^. The copepods were used as the target biological variable 
being by far the most abundant and species-rich groundwater group (—80% by 
individual -based density). 

Meiofauna was extracted by filtering the 20-L samples through a hand net (mesh 
size — 60 ]im). Faunal samples were preserved in an ethyl alcohol solution at 80%. 



Meiofaunal individuals were later sorted, counted, identified to species level, and 
assigned to two ecological categories: obligate (stygobites) or non-obligate ground- 
water species (non-stygobites). Non- stygobites are epigean freshwater species living 
on the springbed, or as inbenthic in sediment interstices (e.g., to avoid predation) or 
are habitat generalists. Some of them are defined as crenobionts^**, i.e. species which 
complete their life cycle in the stable and cold thermal regime of surface spring waters. 

Temperature, pH, specific conductivity at 25 C, and dissolved oxygen (DO) con- 
centration were measured in the field using a multiparametric probe (ECM Multi™; 
Dr. Lange GmbH, Diisseldorf, Germany); calcium concentration was determined by 
titration with EDTA; nitrate, sulphate and phosphate concentrations were measured 
using a HACH DR 2000 spectrophotometer (Hach Co., Loveland, CO, U.S.A.). 
Particulate organic matter (POM) was quantified on 20-L samples, after removal of all 
fauna, oven-drying at 105 C for 24 h, and weighing in anhydrous conditions. The 
dry-weighed samples were ashed at 540'C and re-weighed to determine POM con- 
centration as the difference between dry and ash mass. The grain size composition of 
springbed sediments was measured only once per year at each spring site, as sediment 
transport was negligible throughout the spring system. A more intensive sampling of 
sediments also would have altered habitat features and hence species abundance and 
distribution across sites. The mass of eight granulometric classes was determined at 
each site on the fractional dry-sieving of desiccated samples through a stack of sieves 
according to the Wentworth*^ scale: medium gravel (4-8 mm), fine gravel (2-4 mm), 
very coarse sand (1-2 mm), coarse sand (0.5-1 mm), medium sand (0.25-0.5 mm), 
fine sand (0.125-0.25 mm), very fine sand (0.063-0.125 mm), sUt and clay 
(<0.063 mm). 

Statistical analysis. The sampling design was fully balanced, thus maximising 
statistical power and reliabLhty*^-^". Coleman rarefaction curves run for biological data 
indicated that sampling efforts were adequate for all sampling years (Supplementary 
Table 5, Supplementary Fig. 3). Coleman rarefaction curves were constructed with 
Estimates 9^^ using resampling with replacement based on 50 iterations. 

Between-year comparisons for abiotic and biological variables were tested with 
nonparametric Friedman tests, a block-based variant of the Kruskal-Wallis test, also 
based on ranks, and used when k > 2*^'^^. Sampling sites within sampling years were 
used as blocks, resulting in triplet-based paired comparisons (simultaneous 1997 vs. 
2005 vs. 2012 comparisons). The paired-sample, rank-based approach was chosen 
also to minimise the impact of the inherent between-site variability, especially for 
some variables (e.g., grain size composition of sediments). Sample size (n) was 32 for 
all variables in each year, except for grain size, for which n — S. Friedman tests were 
thus fully reliable for all variables, though the outcome remained slightly liberal for 
the relatively low-replication grain size categories^". However, the Friedman test is 
considered a robust technique as long as the assumption of homoscedasticity is 
j^g^49.5o.52 'j'j^g Friedman test was applied only to those variables with homogeneous 
variances, as tested with a two-tailed Cochran C test^^ {P ^ 0.369 for all abiotic 
variables and major biological groups; Supplementary Table 6). values for Friedman 
tests were determined asymptotically. A Nemenyi paired-sample multiple compar- 
ison was performed for statistically significant Friedman tests to identify pairwise 
year-v5.-year trends. Partly to counterbalance the possible slight liberalism of bor- 
derline significant Friedman and/or Nemenyi results, a Bonferroni correction was 
applied to all tests (whose original significance threshold was set at P < 0.05), leading 
to a Cochran- Friedman -Nemenyi experiment -wide significance threshold of P < 
0.016. Cochran, Friedman, and Nemenyi tests were performed using XLSTAT 
(Addinsoft™ XLSTAT v. 2013.6.01, Paris, France, 2014). 

Data for many individual species of copepods were heteroscedastic (significant 
two-tailed Cochran test: Supplementary Table 6), and could not be tested with the 
Friedman/Nemenyi procedure^"'^^. Year-v5.-year patterns in copepod species distri- 
butions were thus analysed with a test of independence on a contingency table based 
on frequencies of yearly occurrences at the 32 sampling sites (i.e., 32 — the species was 
observed at all sampling sites in that year: Table 1 ). Due to the overall low frequencies 
and patchy species distributions, the original data were coded by a constant X 100 
factor^^. The test of independence was based on two-tailed Pearson distance (X^p) 
with df — (number of years — 1) X (number of species — 1) — 42. Significant 
departures from table-wide independence were followed by an analysis of adjusted 
Pearson residuals to identify observed frequencies for individual species significantly 
higher or lower than expected frequencies^". Significance for adjusted Pearson 
residual analyses was set at the same level as for the preceding Pearson X^ distance, 
standardised for the non-PERMANOVA investigation-wide statistical elaboration 
(i.e., P < 0.016). The same procedure was applied for the two major ecological 
copepod groups (obligate V5. non-obligate groundwater species) in a 3 years X 2 
groups contingency table (Table 1). Contingency tables were analysed with 
Addinsoft™ XLSTAT. 

A permutational multivariate analysis of variance (PERMANOVA)^* was used to 
test for differences in copepod abundances before (1997 and 2005) and after (2012) 
the earthquake. Post-hoc f-tests were used to test for differences between levels within 
the factor. Significance threshold for PERMANOVA analyses was set at P < 0.05. 
Following significant PERMANOVA results, patterns of obligate groundwater 
copepod assemblages before and after the earthquake were examined using non- 
metric multidimensional scaling (nMDS) incorporating the Bray-Curtis similarity 
coefficient. SIMPER analysis^^ was then used to identify the obligate groundwater 
species that mostly accounted for pre- V5. -post-earthquake differences in abundance. 
In both nMDS and PERMANOVA, data were square-root transformed to weight the 
contribution of common and rare species, and a dummy variable of 1 was added to all 
samples to allow the inclusion of otherwise empty samples (zero abundance). A 
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Levene test for homogeneity of variances was carried out before PERMANOVA 
(Supplementary Table 3). All multivariate analyses were done using PRIMER v6 
(PRIMER-E, Plymouth, UK, 2006) and R software {R Development Core Team, R 
Foundation for Statistical Computing, Vienna, Austria, 2008). 
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